function [vnez,svnez,rep] = GNSSlinfit(Mtnez,optinit,optlin,optfig)
%GNSSlinfit - weighted/unweighted least squares linear fit of GNSS time series
%
%   [VNEU,SVNEU,REP] = GNSSlinfit(MTNEU)
%   [VNEU,SVNEU,REP] = GNSSlinfit(TSDATA)
%   [VNEU,SVNEU,REP] = GNSSlinfit([])
%   [VNEU,SVNEU,REP] = GNSSlinfit
%
% Options about the first argument of GNSSlinfit:
%
% - It is a 4-columns coordinate matrix of a GNSS station, MTNEU, where
%       MTNEU = [t NEU] = [t North East Up].
%   In this case the function computes the unweigthed least square linear 
%   fit of each of the time series North, East, Up with the time variable 
%   t. The elements of the 3-by-2 output matrix VNEU are such that
%       VNEU(k,:) = NEU(k,1) + NEU(k,2)*t,   k = 1,2,3.
%   The 2-by-3 output matrix SVNEU provides the standard deviations (STDs)
%   of the computed coefficients, where SVNEU(h,k) is the STD of VNEU(h,k). 
%   In this case, if the data fitting is successful, the output REP is 1.  
%   If the data fitting cannot be performed, VNEU and SVNEU are NaN 3x2 
%   matrices and REP = 0.
%
% - It is a 7-columns matrix
%       MTNEU = [t NEU SNEU] = [t North East Up SNorth SEast SUp],
%   where SNEU = [SNorth SEast SUp] is the matrix of standard deviations
%   corresponding to NEU. 
%   Since in this case the uncertainties are available, also the weighted 
%   linear fit can be performed. The user can interactively choose the 
%   approach to LS fit (weighted or unweighted). If the weighted option is 
%   chosen and the data fitting is successful, it is REP = 2. If the choice 
%   is unweighted, it is REP = 1 as above. 
%
% - It is a tsData object TSDATA. In this case, t, NEU = [North East Up]
%   and SNEU = [SNorth SEast SUp] are taken from TSDATA: t = TSDATA.t, 
%   North = TSDATA.n, East = TSDATA.e, Up = TSDATA.v, SNorth = TSDATA.sn, 
%   SEast = TSDATA.se, SUp = TSDATA.sv. All functionalities are as above.
%
% - It is undefined or empty. In this case, a combo box allows the
%   interactive choice of the MATLAB .mat file with a tsData object. 
%   Please note that two options are possible: the selected file has a 
%   variable only (the tsData object) or has more than one variable and 
%   the field of one of these variables is tsData (case sensitive). 
%   If the data loading is successful, all functionalities are as above.
%
%   [VNEU,SVNEU,REP] = GNSSlinfit(MTNEU,OPTINIT)
%   [VNEU,SVNEU,REP] = GNSSlinfit(TSDATA,OPTINIT)
%   [VNEU,SVNEU,REP] = GNSSlinfit([],OPTINIT)
%
% Options bout the optional input argument OPTINIT:
%
% - If OPTINIT is undefined, empty or false, the function operates as 
%   above.
%
% - If OPTINIT is true and is not a valid GNSS date (see below), an input 
%   dialog box allows the choice of the initial date for the linear fit.
%   in this case, the initial North, East and Up are set to 0 (note that
%   this fact does not imply that the intercept of the straingth line is
%   0).
%
% - If OPTINIT is true and is a valid GNSS date expressed in serial form, 
%   i.e. is not before January 6, 1980 (serial date: 723186)
%
%   [VNEU,SVNEU,REP] = GNSSlinfit(MTNEU,OPTINIT,OPTLIN)
%   [VNEU,SVNEU,REP] = GNSSlinfit(TSDATA,OPTINIT,OPTLIN)
%   [VNEU,SVNEU,REP] = GNSSlinfit([],OPTINIT,OPTLIN)
%
% Options about the optional input argument OPTLIN:
%
% - If OPTLIN is undefined, empty or is neither 1 nor 2, the function 
%   operates and above, i.e. a menu for the choice betweem unweighted or
%   weighted linear fit is shown if the data standard deviations are 
%   available.
%
% - If OPTLIN is 1, the unweighted linear fit is performed, without any 
%   interactive choice;
%
% - If OPTLIN is 2 and standard deviation are available, the weighted 
%   linear fit is directly performed. If OPTLIN is 2 and no data standard 
%   deviations are available, the unweighted linear fit is performed and a 
%   warning message is shown. 
%
%   [VNEU,SVNEU,REP] = GNSSlinfit(MTNEU,OPTLIN,OPTFIG)
%   [VNEU,SVNEU,REP] = GNSSlinfit(TSDATA,OPTLIN,OPTFIG)
%   [VNEU,SVNEU,REP] = GNSSlinfit([],OPTLIN,OPTFIG)
%
% Options about the optional input argument OPTFIG:
%
% - If OPTFIG is undefined, empty of false, no figures are shown.
%
% - If OPTFIG is defined and true, the data and the straight lines fitted 
%   to data are shown.

% G. Teza, 2008, 2021.

if (nargin < 4) || isempty(optfig)
    optfig = 0;
end
optfig = logical(optfig);

if nargin < 3
    optlin = []; 
end
if ~ismember(optlin,[1 2])
    optlin = []; 
end

if nargin < 2 || isempty(optinit) 
    optinit = 0;
end
optinit = logical(optinit);
if nargin < 1 || isempty(Mtnez)
    [filena, pathna] = uigetfile({'*.mat','MAT-files (*.mat)'},...
        'INPUT TSDATA FILE');
    if isequal(filena,0) || isequal(pathna,0)
        [vnez,svnez,rep] = invalidManag('Invalid input file');  % ancillary function
        return
    else
       filenac = fullfile(pathna,filena);
       varfile = load(filenac); 
    end
    fnames = fieldnames(varfile);
    if numel(fnames) == 1
        fname1 = char(fnames);
        Mtnez = varfile.(fname1);
    else
        try
            Mtnez = varfile.tsData;
        catch
            [vnez,svnez,rep] = invalidManag('Invalid input file');  
            return
        end
    end
end

snez = []; 
if isnumeric(Mtnez)         % matrix case
    t = Mtnez(:,1);
    nez = Mtnez(:,2:4);
    if size(Mtnez,2) == 7
        snez = Mtnez(:,5:7);
    end
elseif isa(Mtnez,'tsData')  % object case
    t = Mtnez.t;
    n = Mtnez.n;
    e = Mtnez.e;
    z = Mtnez.v;
    nez = [n e z];
    sn = Mtnez.sn;
    se = Mtnez.se;
    sz = Mtnez.sv;
    snez = [sn se sz];
else
    [vnez,svnez,rep] = invalidManag('Invalid input data');  
    return
end

Inan = isnan(nez(:,1));
nez(Inan,:) = [];
t(Inan) = [];
t1 = t(1);
t2 = t(end);

if optinit
    if (optinit < datenum(1980,1,6)) || (optinit > 2)
        dateinv = datevec(t1); 
        prompt = {'year:','month:','day:'};
        name = 'Initial date for linear fit';
        numlines = 1;
        options.Resize = 'on';
        defaultanswer = {...
            num2str(dateinv(1)),num2str(dateinv(2)),num2str(dateinv(3))};
        answ = inputdlg(prompt,name,numlines,defaultanswer,options);
        ye = str2double(answ{1});
        me = str2double(answ{2});
        de = str2double(answ{3});
        te = datenum(ye,me,de);
    else
        te = optinit;
    end
    Ie = t < te;
    t(Ie) = []; 
    nez(Ie,:) = [];
    t1 = t(1);
    nez(:,1) = nez(:,1)-nez(1,1);
    nez(:,2) = nez(:,2)-nez(1,2);
    nez(:,3) = nez(:,3)-nez(1,3);
else    % to avoid unnecessary calculations in case of unweighted linfit
    Ie = [];
end

n1 = size(nez,1);
if n1 <= 2
    [vnez,svnez,rep] = invalidManag('Unsuccessful linear fit');  
    return 
end
    
if isempty(optlin) && ~isempty(snez)
    optlin = menu('CHOICE OF LINEAR REGRESSION:',...
        'NON-WEIGHTED','WEIGHTED');
elseif isempty(snez)
    if optlin == 2
        warning('No uncertaintes - a non-weighted linear fit is carried out');
    end
    optlin = 1;
end

if optlin == 1

    [vn,sn] = polyfit(t,nez(:,1),1);
    svn = std_fit(sn);  % ancillary function
    vn = flip(vn);      % this because polyfit provides a(2)*t+a(1) 
    svn = flip(svn);
    
    [ve,se] = polyfit(t,nez(:,2),1);
    sve = std_fit(se);
    ve = flip(ve);      
    sve = flip(sve);
    
    [vz,sz] = polyfit(t,nez(:,3),1);
    svz = std_fit(sz);
    vz = flip(vz);      
    svz = flip(svz);
    
    rep = 1;
    
else

    snez(Inan,:) = [];  
    if ~isempty(Ie)
        snez(Ie,:) = [];
    end
    
    [~,an,san,bn,sbn,~,~] = wlinfit(t,nez(:,1),snez(:,1));   
    vn  = [an bn];
    svn = [san sbn];
    
    [~,ae,sae,be,sbe,~,~] = wlinfit(t,nez(:,2),snez(:,2));
    ve  = [ae be];
    sve = [sae sbe];
    
    [~,az,saz,bz,sbz,~,~] = wlinfit(t,nez(:,3),snez(:,3)); 
    vz  = [az bz];
    svz = [saz sbz];
    
    rep = 2;
    
end
vnez  = [vn; ve; vz];     
svnez = [svn; sve; svz];

if optfig
    
    mf = repmat(vnez(:,1),1,2)+vnez(:,2).*repmat([t1 t2],3,1);
    
    figure('units','normalized','outerposition',[0 0 1 1]);
    
    subplot(1,3,1); 
    plot(t,nez(:,1),'.k'); datetick('x'); 
    hold on; plot([t1 t2],mf(1,:),'r','LineWidth',1);
    xlabel('\it t'); ylabel('\it North \rm(m)'); 
    
    subplot(1,3,2); 
    plot(t,nez(:,2),'.k'); datetick('x');
    hold on; plot([t1 t2],mf(2,:),'r','LineWidth',1); 
    xlabel('\it t'); ylabel('\it East \rm(m)'); 
    
    subplot(1,3,3); 
    plot(t,nez(:,3),'.k'); datetick('x');
    hold on; plot([t1 t2],mf(3,:),'r','LineWidth',1);
    xlabel('\it t'); ylabel('\it Up \rm(m)'); 
    
end

%% === ANCILLARY FUNCTIONS ===

function [vnez,svnez,rep] = invalidManag(warningMessage)
% Invalid data management
vnez  = nan(2,3); 
svnez = nan(2,3); 
rep = 0; 
warning(warningMessage); 

function stp = std_fit(s)
% Computes the STDs of the polyfit-based coefficients of the linear fit. 
R  = s.R;  
Rinv = inv(R);
df = s.df;
normr  = s.normr;
covfit = (Rinv*Rinv')*normr^2/df;
stp = [sqrt(covfit(1,1)) sqrt(covfit(2,2))];